Polar Patterns in Active Fluids 



Luca Giomi 1 and M. Cristina Marchetti 2 

1 School of Engineering and Applied Sciences, Harvard University, Cambridge, MA 02138, USA 
2 Physics Department and Syracuse Biomaterials Institute, 
Syracuse University, Syracuse New York, 13^-1130, USA 
(Dated: January 28, 2013) 

We study the spatio-temporal dynamics of a model of polar active fluid in two dimensions. The 
system exhibits a transition from an isotropic to a polarized state as a function of density. The 
uniform polarized state is, however, unstable above a critical value of activity. Upon increasing 
activity, the active fluids displays increasingly complex patterns, including traveling bands, traveling 
vortices and chaotic behavior. The advection arising from the particles self-propulsion and unique 
to polar fluids yields qualitatively new behavior as compared to that obtain in active nematic, with 
traveling-wave structures. We show that the nonlinear hydrodynamic equations can be mapped onto 
a simplified diffusion-reaction-convection model, highlighting the connection between the complex 
dynamics of active system and that of excitable systems. 



I. INTRODUCTION 

Bacterial suspensions, in vitro mixtures of cytoskeletal 
filaments and motor proteins, and migrating epithelial 
cell layers are examples of active fluids composed of in- 
teracting units that consume energy and collectively gen- 
erate motion and mechanical stresses. Due to their elon- 
gated shape, active particles can exhibit orientational or- 
der at high concentration and have been likened to "living 
liquid crystals" pQ. They have been modeled by contin- 
uum equations built by modifying the hydrodynamics of 
liquid crystals to include nonequilibrium terms that ac- 
count for the activity of the system [5HH] , or derived from 
specific microscopic models jTUl E] . 

The theoretical and experimental study of active flu- 
ids has revealed a wealth of emergent phenomena, in- 
cluding spontaneously flowing states [7J [121414] . uncon- 
ventional rheological properties [T51ET] . and new spa- 
tiotemporal patterns not seen in passive complex flu- 
ids (22J 123] • Large-scale swirling patterns and "flock- 
ing" have recently been observed in motility assays con- 
sisting of a highly concentrated suspension of actin fil- 
aments propelled by myosin molecular motor proteins 
tethered to a plane [2U 125] ■ Experiments in bacte- 
rial suspensions [261 EI] and simulations of self-propelled 
hard rods [2"8H3"T] have also shown the formation of po- 
lar high-density clusters traveling trough a low density 
background. In a recent paper one of us and collabo- 
rators examined the spatio-temporal dynamics of active 
fluids with nematic symmetry in two dimensions, demon- 
strating the onset of spatially modulated relaxation os- 
cillations and the close analogy between the dynamics of 
active fluids and that of excitable media 22\ . 

In this paper we consider the spatio-temporal dynam- 
ics of an active suspension with polar symmetry. The 
model is appropriate to describe systems such as bacte- 
rial suspensions, where the particles are swimmers, hence 
head-tail asymmetric or polar. As a result, the active 
fluid can order in a macroscopic polar state, character- 
ized by a nonzero value of a vector order parameter P, 



describing the mean polarization of the system. The or- 
der parameter can be written in terms of a magnitude P 
and a polar director p (with |p| = 1), denoting the di- 
rection of spontaneously broken symmetry in the ordered 
state. While in the nematic state the director n is a head- 
less unit vector because the ordered state is invariant for 
n —> — n, the polar state does not possess this symmetry 
and the polar director p is a true unit vector. As a result, 
the continuum equation for a polar active fluids contain 
terms that are forbidden by symmetry in the case of an 
active nematic. These terms describe swimming of the 
active particles relative to the bulk fluid and yield new 
self-advection contributions to the continuum equations. 

Continuum descriptions of the collective dynamics of 
polar active particles have been discussed before in the 
literature in two cases: for particles moving on a frictional 
substrate ("dry" system) and for particles swimming in 
a fluid (active suspension). The hydrodynamics of the 
dry system is a continuum theory of the Vicsek model of 
flocking, consisting of a collection of self-propelled par- 
ticles moving on a frictional substrate with noisy align- 
ing rules. It was first introduced phenomenologically by 
Toner and Tu [21 [3] and recently derived by explicit coarse 
graining of the microscopic dynamics |32[ 133] - In this 
case the only conserved variable is the density of active 
particles as momentum is not conserved. The theory is 
then formulated entirely in terms of density and polar- 
ization fields and the equations lack Galileian invariance 
due to the presence of the substrate. The hydrodynamic 
equations of active suspensions have been written down 
phenomenologically on the basis of symmetry consider- 
ations and also derived from specific microscopic mod- 
els of a collection of interacting active particles. The 
equations used here were first obtained by Liverpool and 
Marchetti for a suspensions of cytoskeletal filaments with 
crosslinking motor proteins [11] . While in this case the 
filaments may be better described as a "mutually pro- 
pelled", rather than "self-propelled", as the activity re- 
sulting in propulsion arises from the forces exchanged 
among filaments by the active cross-linkers, the form of 
the hydrodynamic equations does not depend on such 
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microscopic details. 

This work focuses on pattern formation in polar active 
suspensions. While the behavior of polar active films in 
a quasi- ID geometry was discussed before in Ref. [14] . 
here we present a systematic study of a two-dimensional 
system with periodic boundary conditions. New results 
include: the derivation of a "phase diagram" for the sys- 
tem that unifies the analysis of instabilities discussed 
previously in the literature in various contexts (Fig. [5]), 
the prediction of the traveling vortices regime shown in 
Fig. |8j the derivation of a low dimensional model that 
captures the onset of the traveling waves regime, and 
the prediction that the transitions between the various 
regimes are discontinuous, allowing for the possibility of 
coexistence and hysteresis. 

We begin by identifying the linear instabilities of the 
ordered state of polar active fluids, summarized in Fig. [2] 
These instabilities have been discussed before in the lit- 
erature, either in the context of dry systems or of active 
suspensions. Here we highlight the connection between 
the behavior of the suspension and the dry limit. To go 
beyond the linear case, we then integrate numerically the 
nonlinear hydrodynamic equations for increasing values 
of activity. In a polar system there exist two classes of 
active terms. The first class, with strength controlled by 
a parameter we call a, is present in both systems of ne- 
matic and polar symmetry and describes active stresses 
induced by coupling of orientation and flow. The sign 
of a identifies active systems that generate contractile 
stresses (a > 0) versus those that generate tensile stresses 
(a < 0). The second class, with strength controlled by 
a parameter we call w, is unique to polar fluids and 
describes a variety of nonequilibrium advective contri- 
butions to the dynamics arising from the particles' self 
propulsion. Both parameters are ultimately controlled by 
the activity of the system, as embodied for instance by 
the rate of adenosine-triphosphate (ATP) consumption 
in motor-filament suspensions or by the forces exerted 
by swimmers on the surrounding fluid in a bacterial sus- 
pension. For this reason, and to reduce the number of 
independent parameters, in most of the following we take 
w = fa, with / a numerical factor of order one, and dis- 
cuss the behavior of the system for increasing values of 
the activity a. Upon increasing a we observe a variety of 
increasingly complex patterns. A regime that is unique 
to polar systems consists of bands extending in the di- 
rection orthogonal to that of mean order (x direction) 
and traveling along x. Although traveling bands have 
been seen in dry systems, either via direct simulations of 
Vicsek models 29J or by numerical solution of the con- 
tinuum theory [321 134j , there are important qualitative 
differences between the traveling bands in dry systems 
and in suspensions. In dry systems the bands appear 
to be solitary waves of ordered regions of high density 
and polarization traveling in a disordered background. In 
suspensions, in contrast, the traveling waves are mainly 
controlled by coupled oscillations in the flow and the di- 
rection of the polarization. They can be understood in 



terms of a simplified reaction-diffusion-advection model 
that highlights the crucial role played by active stresses 
in providing the energy input needed for setting up the 
pattern formation and by advective terms unique to po- 
lar systems in controlling the traveling nature of the pat- 
tern. As the activity is further increased, the traveling 
wave pattern is replaced by more complex spatial struc- 
tures, including oscillating flows with pairs of traveling 
vortices, and eventually chaotic behavior. 

Finally, spatial patterns in active polar fluids in two di- 
mensions have also been discussed by Voituriez et al [35] . 
These authors include a coupling between splay and den- 
sity fluctuations that is present in both equilibrium and 
active polar fluids, but neglect the convective couplings 
unique to active polar systems that are responsible for 
the oscillatory or traveling nature of some of the spatial 
patterns found in the present work. The modulated flow- 
ing phases identified in Ref. [35] are characterized by a 
finite steady flow induced by activity, but do not rep- 
resent traveling patterns. They are "equilibrium-like" 
in the sense that they have a corresponding analogue in 
equilibrium polar systems |36j . When polar convective 
terms are included the modulated patterns are replaced 
by traveling vortices. 

In Sec. [IT] we summarize the hydrodynamic equations 
of a polar active suspension. The homogeneous states of 
these equations and their linear stability are discussed in 
Sec. |III| where we also show the connection between the 
instabilities of the suspension and those of a dry system. 
The results of the linear stability analysis are summarized 
in the phase diagram shown in Fig. [2j In Sec. [iVjwe de- 
scribe the spatio-temporal patterns obtained by numeri- 
cal solution of the nonlinear hydrodynamic equations in a 
planar sample in two dimensions. These include traveling 
bands, traveling vortices and chaotic flow. We also show 
that the nonlinear equations can be mapped onto a sim- 
plified diffusion-reaction-convection model, highlighting 
the connection between the complex dynamics of active 
system and that of excitable systems. We conclude with 
a brief discussion. 



II. HYDRODYNAMICS OF POLAR ACTIVE 
SUSPENSIONS 

We consider a suspension of rod-like active particles 
of length £ and mass M in a fluid. We assume that 
the dynamics of the system can be described by contin- 
uum fields consisting of conserved quantities: the con- 
centration c of active particles, the total density p — 
Psoivcnt + Mc of the suspension, assumed constant, and 
the total momentum density g = pv, with v the flow 
velocity In addition, to describe the possibility of polar 
order we consider the dynamics of the polarization P. 
The concentration and polarization vector of the active 
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particles can be written as [TT] 

c(r,t) = <£*('-*»(*))>, (la) 

n 

P(r,t) = ^y<£ *..(*)*(* -'«(*))) , (lb) 

where r n (t) is the position of the center of mass of the n- 
th swimmer and n (t) is a unit vector along the long axis 
of the swimmer. We consider here uniaxial swimmers, 
which are the simplest type of active polar particles. 

Although the active suspension is a nonequilibrium 
system for which a free energy cannot be defined, it is 
convenient to write the hydrodynamic equations in terms 
of a "free energy" given by 

F = J dr[ia 2 |P| 2 + ia 4 |P| 4 

+ \K8iPjdiPj - w 3 (Sc/c Q )V ■ P] , (2) 

with 5c = c — Co the deviation of the concentration from 
its mean value, cq. The first two terms on the right hand 
side of Eq. ^ determine the homogeneous steady state 
of the system. The parameter a 2 ~ c* — c changes sign 
at a critical concentration c*, while > 0. The third 
term describes the energy cost of spatially inhomoge- 
neous deformations of the order parameters, with K a 
stiffness. For simplicity we have used a one-elastic con- 
stant approximation and equated the elastic constants 
for splay and bend deformations. Finally, the last term 
describes a coupling between concentration and splay de- 
formation that is only allowed in a polar fluid. We have 
neglected polar couplings between splay and fluctuations 
in the magnitude of the order parameter ~ |P| 2 V • P 
that yield advective-type terms and pressure corrections 
in the equation of motion for the order parameter. 

The dynamics of the concentration of active particles 
is governed by a convection-diffusion equation 

3 ( c+V-c(v + Wl P)=V-(DVc), (3) 

where Dij = DgSij + D\PiPj is the anisotropic diffusion 
tensor. The equation for the polarization is given by: 

[d t + (v + w 2 P) • V] Pi + wyPj = XuijPj + j- 1 hi , (4) 

where h = — 5F/5P plays the role of the nematic molec- 
ular field, with: 

h= -(a 2 +a 4 |P| 2 ) P + KV 2 P-w 3 — . , (5) 

c 

Also, 7 is a friction and A a reactive parameter that 
controls the coupling between orientation and flow. As 
in a nematic liquid crystal, order parameter variations 
are coupled to flow gradients, described by the antisym- 
metric and symmetric parts of the rate of strain tensor 
Uij = (diVj + djVi)/2 and = (diVj - djVi)/2. 

Finally, W\ and w 2 are characteristic velocities (taken 
independent of c), proportional to the self-propulsion 



speed of the active particles, and describe self-advection 
of the active particles relative to the fluid. These terms 
are unique to a polar active fluid. The evolution of the 
flow velocity is controlled by the Navier-Stokes equation 

p (d t + v • V) v, = djVij , (6) 

with V • v = 0, as required by incompressibility. The 
stress tensor cry can be written as 

cry = 2r]Uij + a\ 3 + ct?- . (7) 

The first term on the right hand side of Eq. ^ is the 
dissipative part and we have used an isotropic-viscosity 
approximation. The second term is the reversible part, 
given by 

o$ = -mij-^Pihj+Pjh^+^iPihj-PjhiyxSijP-h , 

(8) 

with II the pressure and A another reactive parameter. 
In a nonequilibrium system the coefficients of the various 
terms in the constitutive equation for the stress tensor are 
not in general related to those in the molecular field h, 
as is instead required in equilibrium. Here we have taken 
them to be the same. This simplifying assumption does 
not change the structure of the equation, but only serves 
to limit the number of independent parameters in the 
equations. It is a direct consequence of having written 
the reactive part of the fluxes in terms of the derivative 
of a free energy. Finally, the last term in Eq. ^ is the 
active contribution. It can be written as the sum of two 
terms of different symmetry, er^- = erg + er^ , with 

g% = ac(PiPi + \s i3 P 2 ) , (9a) 
4 = pcidiPj + djPi + 6ijV ■ P) , (9b) 

where erg describes tensile/contractile active stresses 

present in both nematic and polar fluid, while is 
unique to polar fluids and describes active stresses due to 
self-advection and treadmilling [37]. Note that a > cor- 
responds to contractile stresses as generated by pullers or 
by motor-filament systems, while a < describes tensile 
stresses as generated by most swimming bacteria, such 
as E. coli [3S]. In the following we will ignore the polar 
contribution to the active stress and simply incorporate 
the effect of polarity via the advective terms proportional 
to Wi. We will further assume w± = w 2 = 2^3/7 = w. 
We make the system dimensionless by scaling all lengths 
with the length £ of the active particles, scaling time with 
the relaxation time of the director r p — £j/K and scaling 
stress using the elastic stress a = Kjl 2 . 

III. HOMOGENEOUS STEADY STATES AND 
THEIR STABILITY 

There are two homogeneous, steady state solutions of 
the hydrodynamic equations, both with constant density 
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Co. For cq < c* , or 020 = 02(20) > 0, the homogeneous 
solution is the isotropic state, with P = 0. For c > c* , 
or a 2 o < 0, rotational symmetry is spontaneously broken 
and there is an ordered polarized state, with |P| = Pq = 
\J — «2o/ a 40- In the numerical work described below we 
have used the following form for the coefficients a 2 and 
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a-2 
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This choice ensures that Pq ~ ^/c — c* for Co — > c* + and 
P ~ 1 for c > c*. 

Next, we examine the linear stability of the homo- 
geneous states against spatially inhomogcncous fluctu- 
ations in the continuum fields. The isotropic state is al- 
ways linearly stable and will not be discussed further. In 
order to ensure the incompressibility condition, it is con- 
venient to rewrite the Navier-Stokes equation in terms of 
vorticity and stream function tp by expressing: 



dyip 



-d x ip 



The incompressibility condition is then automatically 
satisfied. The vorticity field is given by 



u> = 2u xy = d x v y - d y v x , 



(10) 



and the stream-function if> is related to the vorticity w 
through a Poisson equation: 



(11) 



The Navier-Stokes equation can be written in terms of w 
and ip as 



3 t UJ = TjS7 2 UJ + d 2 x Ty X + d X yTyy ~ 8y X T XX 



y'xy 



, (12) 



where we have separated out the viscous part of the stress 



tensor and defined Tjj as 



nj = cr'-j + afj - -Sij (a 7 kk + a kk ) 



(13) 



To analyze the stability of the polarized state, we 
choose the x axis along the direction of order and 
compactify the notation by introducing a vector ip = 
{c, P x , P y ,u)}. The homogeneous polarized state is char- 
acterized by ip(°' = {co, Pq, 0, 0}. We then introduce the 
fluctuations from the homogeneous values by letting 



¥>(x,t) = ¥>W+e V W(x,t) 

and consider the dynamics of fluctuations for e <C 1. 
To enforce periodic boundary conditions on a square do- 
main, we look for solutions of the form 



v (1 W)= £ 
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where q n — 2im/L and q m = 2nm/L. The Fourier com- 
ponents of the stream-function are related to those of the 
vorticity by: 



nm — 9 
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FIG. 1: Real part of the eigenvalues s± of the matrix aio 
displaying the coupled instability of (<5c, SP X ) for w = 0.5 and 
cq = 0.01. 



With this choice the linearized hydrodynamic equations 
reduce to a set of coupled linear ordinary differential 
equations for the Fourier modes y nm , given by 



dtfnm = A nm ip ri 



(14) 



The matrix A nm can be expressed in the following block- 
form 



(15) 



The explicit expression for the matrix A nm is given in 
Appendix[S] In order for the homogeneous state ipo to be 
stable, the real part of the eigenvalues of the matrix A nm 
must be negative. It is useful to examine the behavior of 
the eigenvalues for (n,m) modes (1,0) and (0, 1), corre- 
sponding to wavevectors longitudinal and transverse to 
the direction of mean order, respectively. In both cases 
(5c, P x ) fluctuations decouple from (P y ,w) fluctuations 
and the dynamics can be examined by solving two inde- 
pendent 2x2 eigenvalue problems. 

We first consider the (1, 0) or longitudinal modes. The 
matrix 610 = 0, hence Aio is block diagonal. The cou- 
pled dynamics of c and P x fluctuations (note that 5P X 
describes fluctuations in the magnitude of the order pa- 
rameter) is controlled by the eigenvalues of the 2x2 
matrix aio, which does not depend on the active param- 
eter a, but only on the advective velocity w. The (10) 
eigenvalues are given by 



with: 



--(C + 2iqwP ) 



2w 2 q 2 — SiqwPQ 



c c* 



co 



(16) 



C = 2\c Q -c*\ + (K + D)q 2 , 

q = 2-k/L and Dq = P>\ = D. The fluctuations of con- 
centration and order parameter magnitude in the longitu- 
dinal (1,0) exhibit unstable growth for arbitrarily small 
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FIG. 2: (color online) Phase-diagram in the plane (co — c* , a) 
for w = a and q = 0.2, A = 0.5, t] = D = £>i = 1. The 
labeled regions correspond to the homogeneous steady state 
(HSS), the spontaneous flow regime (SF), the traveling waves 
regime (TW) and the region of longitudinal traveling waves 
of Vicsek-type ( V) . 



values of w as the mean field order-disorder transition is 
approached from above, i.e., cq — > c* + . The real part of 
the eigenvalues of the matrix a%o is shown in Fig. [T] This 
instability does not couple to the flow and is precisely the 
instability obtained in Refs. [32j [34] for dry systems con- 
sisting of self-propelled particles with aligning interac- 
tions moving on a frictional substrate. It occurs for both 
tensile and contractile systems, but it is unique to po- 
lar active fluids. The instability only occurs in a narrow 
range of densities above the putative mean-field contin- 
uous transition point c* between polarized and isotropic 
states, as shown in the stability phase diagram of Fig. [2] 
It has been argued that it may be associated with the 
renormalization of the order of the transition, which has 
been shown to be discontinuous in recent simulations of 
Vicsek-type models 29 . We note that the magnitude 
of the order parameter P x is not a hydrodynamic vari- 
able in the strict sense as its decay rate has the finite 
value — 2|a2o| at large wavelengths. At the transition, 
however, 0,20 = and fluctuations in P x become long- 
lived, driving the instability. In the absence of noise, the 
numerical solution of continuum equations obtained by 
coarse-graining the Vicsek model has yielded in this re- 
gion of the parameters a spatially inhomogeneous state 
consisting of solitary waves or bands of ordered regions 
traveling in a disordered background [32] ■ We refer to the 
region of parameters where this linear instability occurs 
as longitudinal traveling Vicsek-type waves (V) regime. 

The coupled dynamics of fluctuations in the transverse 
component P y of the order parameter, describing director 
fluctuations, and vorticity w yields additional instabili- 
ties, both in the (1, 0) and (0, 1) directions, controlled by 
the d nm matrix. The first unstable modes are the longi- 
tudinal mode (1, 0) and the transverse mode (0, 1). These 
modes becomes unstable for negative and positive values 
of the combination a(l — A), respectively. These instabil- 



ities arise from the coupling of director deformations and 
flow and occur even for w — 0, corresponding to systems 
with nematic symmetry. For a(l — A) > 0, transverse 
modes (0, 1) go unstable, corresponding to coupled splay 
deformations of the director and fluctuations of the v y ve- 
locity. For A < 1, the instability occurs above a positive 
critical value of a given by 



q 2 [^ + P 2 (l-\f 
2c P 2 (l- A) 



(17) 



More precisely this instability occurs for tumbling (|A| < 
1) or disk-like (A < —1) pullers for a > aoi > and for 
rod-like (A > 1) pushers for a < aoi < 0. The interplay 
between the parameter A controlling the flow alignment 
properties of an orientationally ordered suspension and 
the activity a in controlling the instabilities of active flu- 
ids was discussed in Ref. [H] and we refer to that work 
for details. The critical value aoi of transverse or bend 
instability does not depend on w and is identical to the 
value obtained for active nematic, with the replacement 
Pq — > S, where S is the magnitude of the nematic or- 
der parameter. This instability has been obtained before 
in the literature [3] and is commonly referred to as the 
generic instability. For a(l + A) < 0, corresponding to 
tensile active stresses as generated by pushers, such as E. 
coli, longitudinal modes (1,0) go unstable, corresponding 
to coupled bend deformations of the director and fluctua- 
tions in the v x velocity. This instability occurs for A > — 1 
and a < aio < and for A < — 1 and a > aio > 0, with 



«io 



4 m 2 (l + rj) 2 + Pg[Arjw 2 + q 2 (l + rj) 2 (l + A) 2 



2 Co P 2 (l + A)(l + ?? ) 2 



(18) 



Notice that the critical value aio depends on w (a finite 
w in this case actually suppresses the instability). When 
w = 0, this coincides with the generic instability dis- 
cussed in [3] for pushers and active nematic. For finite w, 
however, the modes that goes unstable is a propagating 
wave, suggesting that the system may support traveling 
waves solutions in this region of parameter. This region 
is labelled at TW'm Fig. [2] Finally, the splay/bend insta- 
bilities obtained here are simply the 2D generalization of 
the instability to a spontaneously flowing state discussed 
in Refs. [T^] and [2] for nematic and polar active fluids, 
respectively, in a quasi-one-dimensional strip geometry. 
Their origin lies in the long-range nature of the hydrody- 
namic interactions between active swimming particles, as 
shown in Ref. [39 . Some additional detail on the eigen- 
values of the d nm matrix can be found in Appendix [B| 



A. Relation to Dry Systems 

For completeness we compare the behavior of the ac- 
tive suspension considered here to that of active particles 
on a frictional substrate, referred to so far as dry systems. 
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In this case the only conserved field is the concentra- 
tion of active particles and the continuum equations are 
written solely in terms of concentration and polarization 
fields. The momentum of the system is not conserved and 
there is no equation for the flow velocity 40J. The low 
density longitudinal instability of the ordered state ob- 
tained near the mean-field order order-disorder transition 



is unique to polar active fluids with aligning 



interactions, corresponding for instance to coarse-grained 
versions of the Vicsek model, and does not couple to the 
flow velocity. It is present in both dry systems and sus- 
pensions. It was discussed for the case polar particles on 
a substrate in Refs. [32] and [53] , 

In contrast, the splay and bend instabilities obtained 
at high density are due curvature-induced currents aris- 
ing from the coupling of director deformation and flow. 
One can obtain the limit of dry systems from our gen- 
eral equations for a suspension by adding a drag — £v to 
the right hand side of Eq. ([6]), describing the coupling 
frictional coupling to a substrate. In the limit of large 
friction one can then neglect inertial terms and write 
the solution of Eq. ^ to lowest order in the gradients as 



-8,11 



(19) 



where the active stress <r"- is given in Eq. (9a) and we 



must satisfy V • v = 0. If Eq. ( 19 1 is used to eliminate 



v from the equations for concentration and polarization, 
one then recovers both the splay and bend instabilities. 
Note that when Eq. ( 19 ) is used to eliminate the flow 



velocity in favor of polarization and density the terms re- 
sponsible for the splay and bend instabilities are of order 
(V 2 P 3 ). These terms have been neglected in the contin- 
uum models of polar fluids on a substrate discussed in the 
literature, which explain why such instabilities were not 
obtained to linear order [IT]. When considering the full 
nonlinear equations similar effects are, however, provided 
by other nonlinearities. 



IV. NONLINEAR SPATIO-TEMPORAL 
PATTERNS 

In this section we report some results obtained from 
a numerical solution of Eqs. ([3]), Q and ([6]). To ex- 
plore different regimes at increasing level of activity we 
set w — 0.1a and use a as a variable parameter. Eqs. ([3]), 
Q and Q are then integrated using a vorticity/stream- 
function finite difference scheme on a collocated grid of 
lattice spacing Ace = Ay = 0.078. The time integration is 
performed via a fourth order Runge-Kutta met hod with 



time step At = 10 . As mentioned in Sec. 
vorticity /stream- function method requires one to solve a 
Poisson equation at each time step in order to calculate 
the components of the flow velocity. This was efficiently 
preformed with a U-cycle multigrid algorithm. As initial 
configurations we consider a homogeneous system whose 
director field was aligned along the x axis subject to a 
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FIG. 3: (color online) The velocity field (left) and the po- 
larization direction (right) superimposed to a density plot of 
the concentration and the polarization magnitude for a = 1.5 
(top) and q = —2.1 (bottom). The colors indicate regions of 
large (red) and small (green) density and large (red) and small 
(blue) polarization. For positive values of a the homogeneous 
state undergoes a splay instability that leads to the formation 
of two stationary bands. For negative a a bending instability 
results in the appearance of vertical bands traveling at speed 
vq ~ w ~ a (hence from right to left). 



small random perturbation in density and orientation. 
Thus c = c + e, 9 = e, P = y/(c- c*)/(c + c*) and 
v x = v y = 0, where e is a random number of zero mean 
and variance (e 2 ) = 10~ 2 . The equations where then in- 
tegrated from t = to t = 10 3 , corresponding to 10 6 
time steps. For all the numerical calculations described 
in this section we used rj = c* = Dq = D% = 1, A = 0.1 
and L = 10, in the units described in Sec. [TTJ . Upon 
varying the activity parameter a five different regimes 
have been observed: 1) homogeneous stationary state; 
2) spontaneous steady-flow; 3) oscillating flow with or- 
thogonal traveling bands; 4) oscillating flow with pairs of 
traveling vortices and 5) chaotic flow. These regimes are 
described below in more detail. 



Homogeneous State and Spontaneous Steady 
Flow 



For small values of |a|, the system quickly relaxes to a 
homogeneous state with the polarization vector aligned 
along the x direction. The concentration is equal to the 
initial concentration cq and the polar order parameter 
equates its equilibrium value P a — \J (c — c*)/(co + e*). 
Upon raising a above the critical value (17) , the homo- 
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FIG. 4: (color online) The quantities c, v y , P x and P y along 
the x axis in the traveling waves regime for a = —2.1. The 
yellow region indicate the band visible in the bottom panel of 
Fig. [3] The slant of the yellow region indicate the direction 
of motion of the band. 
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FIG. 5: (color online) The quantities c, u), P x and P y at the 
point x — y = L /2 as a function of time in the traveling waves 
regime for a — —2.1. 



geneous state becomes unstable to a non-homogeneous 
flowing state (top of Fig. [3|. The structure of this state 
is the same of active nematic fluids and consists of two 
bands flowing in opposite directions [55]. The direction 
of the streamlines is dictated by the initial conditions 
which, in this case, favor a flow in the x direction. This 
solution has been intensively discussed in the literature 
of active matter (see for instance [T2T - TH] ) and will not be 
commented here any further. 



B. Traveling Waves 



flow, thus forming traveling bands across the x direc- 
tion. The polarization direction p, on the other hand, 
oscillate while keeping the mean orientation orthogonal 
to the bands. 

The physical origin of the traveling waves in this model 
of active polar fluids relies on the fact that the polar- 
ization field is simultaneously a broken-symmetry vari- 
able, a mechanical stress and a physical velocity. Being 
a broken-symmetry variable, the polarization relaxes to- 
ward a non-zero value in the characteristic time scale r p 
(see Sec. [nj). Correspondingly, shear stress a® y ~ aP x P y 
is injected in the system at rate t" 1 ~ a. When r a rs r p 
this residual stress is accommodated in the passive struc- 
tures leading to an elastic distortion and flow as discussed 
in Sec. |III| Unlike nematic fluids, here the existence of a 
direction of broken-symmetry results in a directed motion 
of particles with velocity v a ~ wP, hence the pattern 
produced by the interplay between active stress, liquid 
crystalline elasticity and flow is advected along P result- 
ing in the traveling waves described above. Note that the 
pattern generated by the splay instability (top of Fig. [3 1 , 
is not accompanied by the formation of traveling bands 
since the system is translational invariant along the po- 
larization direction, hence there is no advective trans- 
port: i.e. P V = P x d x = 0. 

To elucidate the properties of this solution, it is use- 
ful to reduce the hydrodynamic equations to a tractable 
form. In the following we will see that under drastic 
but reasonable simplifications, the hydrodynamic equa- 
tions presented in Sec. [TT] map into a reaction-diffusion- 
advection system that admits traveling waves solutions. 
Diffusion-reaction systems have played a fundamental 
role in the study of biological waves and pattern forma- 
tion, especially in the context of morphogenesis since the 
seminal 1952 paper by Turing [3SJ. 

As a starting point we notice that at the onset of the 
traveling bands regime, the variations in c and P x are 
only slight with respect to those in oj and P y (see Figs. 
[3] and [5]) . This suggests that the dynamics of the travel- 
ing bands is dominated by the latter two fields, as also 
indicated by the linear stability analysis. Thus we ap- 
proximate c and P x as constants and focus our attention 
on the equations for uj and P y : 

[8 t + (V + VjP) ■ V]Py = X(Uy X P X + UyyPy) ~ LO yX P X + fly . 

Next, we notice that v x = and that there is no de- 
pendence on the y coordinate. Incompressibility ensures 
then: u TT — u m — 0. Morover u rv — u VT = lu 



-u> yx — uj/2. The previous expression and Eq. (12 1 sim- 



plify then to: 



For negative values of a larger in magnitude than the 
critical value (18) (a s» —1.19 for our choice of param- 
eters) traveling waves form in the system. As shown in 
Fig. [3] the waves consist of a shear-flow along the y direc- 
tion, that translates at constant speed along x: v x = 
and v y = v y (x — v t) with v the wave velocity. The 
concentration c and the polarization field P follows the 



d t Py = -P X (1 + A)W - WP X 8 X Py 



[a-2 



ai (P x + Py)]P y + dlP, 



y ' 



d t oj = rjdlu + d 2 x T yx 



(20a) 
(20b) 



Without loss of generality, we take 04 = 1 and + 
claP x = —a- Now, in order to understand the origin of 
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onset of the transition, we can construct an approximate 
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FIG. 6: (color online) The frequency v of the waves and the 
magnitude of P y as a function of a in the traveling waves 
regime. 



the traveling bands, one does not need to consider the 
full stress tensor T yx , but only the active contribution 
T yx — aP x P y , which is responsible for injecting shear 
stress inside the system at the rate a. Thus, taking for 
simplicity w — a and transforming a — ¥ aP x we arrive 
to the following set of equations: 



P + aP 1 = {a- P 2 )P + P" 
U) = rjoj" + aP" , 



(21a) 
(21b) 



where the dot and the prime indicate respectively the 
time and space derivative, b = P x (l + A)/2 and we have 
renamed P = P y for shortness. Eqs. (21) have the clas- 
sical form of a diffusion-reaction-advection system. The 
cubic nonlinearity is the trademark of excitable dynami- 
cal systems and, as discussed in Ref. [52] , gives the active 
fluid the behavior of an excitable medium. The advec- 
tion term aP', on the other hand, is an exclusive feature 
of active polar systems and is ultimately responsible for 
the existence of traveling waves. 

For small values of a the solution of Eqs. (21 ) is the 
homogeneous state (P, uj) = {±y/a, 0). A small pertur- 
bation of the homogeneous state of the form SP(t) = 
Pk{t)e lkx and 5oj(t) = oj k (t)e lkx evolves according to the 
equations: 



± ( SP k 
dt \ Suj k 



—2a — k 2 — ika 
-ak 2 



b 

-rjk 2 



SP k 
5u k 



Thus, upon lowering a below the critical value a c ~ 
~2ar\jb (for small k) the homogeneous state becomes un- 
stable and the solution consists of a wave train traveling 
along the negative x direction at constant velocity. In or- 
der to describe the behavior of the traveling waves at the 



solution of Eqs. (21 ) using the method of harmonic bal- 



ance. We then look for a solution of (21 ) of the following 
form: 

P(x,t) = P k {t)e ikx + P k (t)e- ikx , 

u(x,t) = u k {t) lkx + uj k (t)e~ lkx , 

with the bar indicating the complex conjugation. Sub- 
stituting this in (21) and taking only the lowest wave 



numbers yields the following system of ordinary differen- 
tial equations for the amplitudes: 



\k 2 - a + ika)P k - iP k P k 



buj k 



Pk = - 

Pk = -(k 2 -a- ika)P k - 3P k P k 

w fe = -k 2 (i]Lj k + aP k ) , 

ui k = -k 2 (T]Uj k + aP k ) . 

These equations admits solutions of the form: 

Pk{t) = \P k \ e^ t+M , u k {t) = \ Uk \ e^ t+ ^ . 

Without loss of generality we can take <fri — and 
4>2 = 4>. Then, substituting this form in (22) we find, 



(22a) 
(22b) 
(22c) 
(22d) 



after some manipulations, the following set of algebraic 
equations for the quantities is, 0, \P k \ and \uj k \: 



• 




3|P fc | 3 + {k 2 -a)\P k \-b\u k \cos0 
ak 2 \P k \ + {-qk 2 cos0 — v suv(j))\u) k \ 
(ka + v)\P k \ - b\uj k \ sin0 = , 
v cos <j) + r]k 2 sin (f> = . 
Solving this equations to leading order in a gives: 

/m . - — . (24) 



(23a) 
(23b) 
(23c) 
(23d) 



a 
krj 



\Pk\ 



I rj(a — k 2 ) — ba 



'Srj 



' r](a — k 2 ) — ba 
3r] 3 



The larger spatial mode k — 2tt/L thus evolves in time 
through traveling waves of the form: 

P(x,t) w 2|P fc |cosfc(a; - at) 

uj(x, t) w 2\uj k \ cos k(x — at — a/r/k 2 ) 

Fig. [7] shows a plot of the amplitudes \P k \ and \uj k \ for 
the largest wavelength corresponding to k = 2t:/L in 
Eqs. p4| together with those obtained from a numerical 
solution of the simplified equations (21 ). 



The combined linear stability and harmonic balance 
analysis show that the Hopf bifurcation that leads to the 
formation of traveling waves is subcritical with bistable 
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FIG. 7: (color online) The amplitudes of the traveling waves 
resu lting from the numerical solution of the simplified model 
(21 1 together with the predictions of the harmonic balance 
analysis (24 1 for the largest spatial mode k = 2-n/L (solid 
lines). At a « —2ar\jb the system undergoes a subcriti- 
cal Hopf bifurcation with a bistability region in the range 
—2ar\jb < a < 0. The red squares and blue/green circles cor- 
respond respectively to the stable homogeneous steady state 
and the stable large-amplitude limit cycle. 



region in the range —2ar\jb < a < 0. As the activ- 
ity parameter a is lowered below zero, the system goes 
through a subcritical bifurcation where there is a stable 
stationary state and a stable large-amplitude limit cycle 
with an unstable periodic orbit acting as the separatrix 
between the two states. The existence of bistability at 
small negative values of a could have significant effects 
in the collective motion of active polar suspensions such 
as the presence of hysteresis when the amount of activity 
is cycled across the bifurcation point. 



C. Traveling Vortices and Chaos 



FIG. 8: (color online) The velocity field (left) and the polar- 
ization direction (right) superimposed to a density plot of the 
concentration and the polarization magnitude at two different 
times for a = — 5. The flow is characterized by two vortices 
and two stagnation saddle points traveling across the system 
from left to right. 
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FIG. 9: (color online) The quantities c, u, P x and P y at the 
point x = y = L/2 as a function of time in the traveling 
vortices for a = —5. 



If we further decrease a, the traveling waves described 
in the previous section become unstable. In the new 
regime, the concentration is distributed in the form of 
kidney-shaped clusters that travels along the positive x 
direction at constant speed (Fig. fright). The flow field, 
in turn, consists of a pair of oppositely spinning vortices 
whose centers also move at constant speed following the 
concentration clusters as well as two saddle stagnation 
points due to the toroidal topology of the domain (Fig. 
H left). The order parameter appears concentrated along 
"valleys" whose curved conformation partially resembles 
the banded structure of the previous solution. Once again 



the relative variation in polarization is larger than that 
in concentration so that the value of the order parameter 
at the bottom of the valleys is almost 50% lower than 
along the ridges. 

The formation of the vortices is associated with a 
break-down of the translational invariance along the lon- 
gitudinal direction of the bands. At the onset of vortex 
formation the vorticity function for the spatial mode of 
largest wave-length is given heuristically by: 

Wfe(x, t) ~ lo±_ cos k(x — aP x t) — wy cos ky (25) 
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whose corresponding velocity field consists of two vortices 
and two saddles traveling along the x direction with ve- 
locity aP x . From the numerical data shown in Fig. [8] 
we see that P x is everywhere negative, hence the vortices 
move toward the positive x direction. In contrast with 
the traveling waves, where the polarization vector oscil- 
late around its initial orientation (x in this case), the 
formation of traveling vortices is anticipated by an inver- 
sion of the polarization vector (i.e. p ~ — x). The origin 
of this inversion is likely due to the short time scale at 
which the active stress is initially injected in the system, 
which makes the polarization field undergo dramatic dis- 
tortions before it is able to settle on a limit cycle. 

Upon further decreasing a, the traveling vortices be- 
come unstable and the system undergoes a transition to 
a chaotic regime. The polarization valleys now continu- 
ously form and disappear in unpredictable fashion. The 
flow still exhibits larger vortices, which however move 
chaotically across the sample with variable direction and 
magnitude. Interestingly, the polarization direction is 
continuously distorted and does not posses regions of 
partial alignment. This differs from the spatio-temporal 
chaos found active nematics [25] where the chaotic flow 
is organized in grains of uniform orientation separated by 
chaotically moving grain boundaries. The nonlinear spa- 
tiotemporal patterns described in this section are sum- 
marized in the phase-diagram of Fig. 11 The Vicsek- 



type traveling waves shown in the linear phase diagram 
of Fig. [2] are also obtained from the numerical solution 
of the nonlinear equations, but at a much lower density 



that shown in Fig. 11 



V. DISCUSSION 

We have analyzed the complex spatio-temporal pat- 
terns of an active polar fluid in two dimensions. The 
continuum model considered is appropriate to describe 
systems such as bacterial suspensions, consisting of self- 
propelled particles in a solvent. The self-propulsion of 
the active units yields a number of advective terms in the 
hydrodynamic equations that are unique to polar fluids. 
These advective terms are responsible for the traveling 
nature of the large scale structures, which include travel- 
ing stripes oriented normal to the direction of propaga- 
tion and traveling vortices. These structures are remark- 
ably stable in a wide region of parameters and resemble 
the patterns seen recently in actin motility assays at high 
actin density, suggesting that systems may indeed have 
an underlying polar symmetry. We stress that the tran- 
sitions between regimes characterized by different stable 
large-scale structures can only be obtained by allowing 
for spatial variations of the order parameter. A novel 
feature of the present work is that the transitions be- 
tween these various regimes appear to be discontinuous, 
allowing for the possibility of coexistence and hysteresis. 
Discontinuous transitions between disordered and and or- 
dered states have also been seen in dry active systems via 
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FIG. 10: (color online) The velocity field (left) and the po- 
larization direction (right) superimposed to a density plot of 
the concentration and the polarization magnitude for a — —6. 
The flow consists of large vortices that move chaotically across 
the system with variable direction and magnitude. 
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FIG. 11: (color online) Phase diagram of the nonlinear spa- 
tiotemporal patterns: homogeneous steady state (HSS), trav- 
eling waves (TW), traveling vortices (TV) and chaos (CH). 
The dots have been obtained from the numerical integration 
of Eqs. (3), @ and {B) for A = 0.1, r) = D = Dt = 1, 
L — 10 and w — 0.1a. The solid blue line separating the 
homogeneous steady state from the traveling waves regime is 
given by Eq. (18 1. The dashed line separating the traveling 



waves and the traveling vortices regime has been inferred from 
the data. The color gradient between the traveling vortices 
and the chaotic regime indicates a fuzzy boundary line. 
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numerical studies of the Vicsek model [55] . Understand- 
ing the nature and order of the transition in both sus- 
pensions and active particles on substrates will, however, 
require further investigation. 

The present analysis is limited to two dimensional sys- 
tems and may be applicable to situations like motility 
assays where the system is quasi-2D. One should, how- 
ever, be cautious about generalizing the results found in 
2D to 3D systems, as new patterns may arise one the po- 
larization field is allowed to rotate out of the plane. Very 
recent work by Fielding et al. [23] has indeed shown that 
the behavior of a sheared active nematic in two dimen- 
sions can be quite different from that in one dimension. 

Our work, together with recent work on pattern for- 
mation inactive nematic, suggests a strong analogy be- 
tween the spatio-temporal dynamics of active fluids and 
the behavior of excitable systems. In both cases the 
complex behavior obtained from the nonlinear hydrody- 
namic equations can be reproduced qualitatively by sim- 
plified lower dimensional models that resemble reaction- 
diffusion or reaction-diffusion-advection equations famil- 



iar from the study of excitable systems. This mapping 
may provide a useful framework for understanding and 
classifying active systems. 



Acknowledgments 

MCM was supported by the NSF on grants DMR- 
0806511 and DMR-1004789. LG was supported by the 
Harvard-NSF MRSEC, the Harvard-Kavli Nano-Bio Sci- 
ence and Technology Center and the Wyss Institute. We 
thank Aparna Baskaran and Paolo Paoletti for illuminat- 
ing discussions. 



Appendix A: Appendix A 

In this appendix we give the explicit expressions for 
the elements of the matrix A nm given in Eq. (15 1. The 
various block elements are given by 



-iwq n c 

-(a> 2 + a' 4 P 2 )P -K (ql + ql) - twq n P Q ~ 2\a 20 \ 



(Al) 



-iwq m co 

ri q m q n XP 



2 c 



L i 



21 f~i22 

nrn ^ 'ram 



with: 



r- 



aPo (<Zm 



q n Qr, 



2<7m9nAc* Pq 

c + c* 



iq m wP a (l - A) (ql + ql) 
4c 



Cram = 2c o P o [q m q n ^ + a (ql + q n q m - ql)] + q m q n XP [(ql + ql) - 2c* 



and: 



^rj.m. 



- [0?m + ql) + knwP ] P °^ (1 2 (q£+qi) X+1) ^ 

CoaPo (ql - ql) - f (q 2 m + q 2 n ) [(1 - \)q 2 m + ql(X + 1)] -7? (ql + q 2 n ) 



where 020 = 02 (co) = c* — cq < 0, Do — D>\ — D and the prime denotes a derivative with respect to concentration, 
i.e., a' 2i = {da 2 A/dc) c=Cg , which gives -(a' 2 + ci^Pq) = 2c* /(c + c*). 



Appendix B: Appendix B 

Analytical expression for the values of a associated 
with these instabilities can be found by diagonalizing 



the matrix d nm . For both the longitudinal (10) and the 
transverse (01) directions, it is the coupled modes of P y 
and oj that go unstable. For a > the modes go un- 
stable along the 01 direction in this case P y corresponds 
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to splay deformations and the vorticity describes fluctua- and: 
tions in the v x velocity. For a < the modes go unstable 
along the 10 direction in this case P y corresponds to bend 

deformations and the vorticity describes fluctuations in / — q 2 |(1 — A)Po 

the v v velocity. Both these instabilities exist also when ^oi — I , n ri 2ci w i 2 

w = 0. lhe matrices dio and doi are given by: x A 

-q 2 -iqwP §(1 + A)flj 



dio = [ ) with q = 2n/L. The corresponding eigenvalues are: 

'"[ 2 < 



-(7 2 Po[k 2 (l + A)+ac ] -rjq 2 



Aio = - ig |?(»7 + 1) + *wP ± ^(77 - l) 2 - P 2 K + q 2 (l + A) 2 + qco(1 + A)] - 2iqwP {v - 1) 
A i = -\q Uv + 1) ± - l) 2 + P 2 (l - A)[2ac - <z 2 (i _ A )]| 
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